home *** CD-ROM | disk | FTP | other *** search
/ Developer CD Series 1999 June: Reference Library / Dev.CD Jun 99 RL Disk 1.toast / Technical Documentation / Develop / develop Issue 28 / develop Issue 28 code / CurveLayout / AccurateShapesLibrary / ComputeLength.c < prev    next >
Encoding:
Text File  |  1996-09-22  |  9.7 KB  |  328 lines  |  [TEXT/CWIE]
  1. //
  2. //
  3. //         File:    ComputeLength.c
  4. //
  5. //        Joseph Maurer, April 1996
  6. //
  7. //
  8.  
  9.  
  10.  
  11. #include <stdio.h>
  12. #include <Types.h>
  13. #include <fp.h>
  14.  
  15. #include "ComputeLength.h"
  16.  
  17.         /* better than 1/2 pixel */
  18. #define kEpsilon        0.1
  19.  
  20. //-----------------------------------------------------------------
  21. struct CurveCoefficients {
  22.     double a0, a1, a2;    // x(t) = a0 + t * (a1 + t * a2);
  23.     double b0, b1, b2;    // y(t) = b0 + t * (b1 + t * b2);
  24.     Boolean isLinear;    // true if all three points are colinear.
  25. };
  26. typedef struct CurveCoefficients CurveCoefficients, *CurveCoeffPtr;
  27.  
  28.  
  29. //------------------------------------------------------------------
  30. // A (16.16) Fixed number is a 32-bit long integer "n" representing
  31. // a decimal value of (n / 65536).
  32. #define kFloatToFixed    65536.0
  33.  
  34. #define MULTIPLY_DOUBLE_WITH_FIXED(d, n)    (Fixed)((d) * (double)(n))
  35. //------------------------------------------------------------------
  36. void CurveToLine(gxCurve *theCurve, gxLine *theLine);
  37. void CurveToLine(gxCurve *theCurve, gxLine *theLine)
  38. {
  39.     // function takes a curve that is really a line and makes a line structure out of
  40.     // it.  It checks the 3 combinations of the three points to find the two furthest apart
  41.     //    because it is not necessarily the first and last points.
  42.     //    use doubles to avoid overflow.  Function assumes the curve was all colinear points.
  43.     double        x1, y1, x2, y2, x3, y3;
  44.     double        d1, d2, d3;
  45.     
  46.     x1 = theCurve->first.x / kFloatToFixed;
  47.     y1 = theCurve->first.y / kFloatToFixed;
  48.     x2 = theCurve->control.x / kFloatToFixed;
  49.     y2 = theCurve->control.y / kFloatToFixed;
  50.     x3 = theCurve->last.x / kFloatToFixed;
  51.     y3 = theCurve->last.y / kFloatToFixed;
  52.     
  53.     d1 = (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1);            // line: first-control
  54.     d2 = (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1);            // line: first-last
  55.     d3 = (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2);            // line: control-last
  56.     
  57.     if ( (d1 > d2) && (d1 > d3) ) {        // d1 longest
  58.     
  59.         theLine->first.x = theCurve->first.x;
  60.         theLine->first.y = theCurve->first.y;
  61.         theLine->last.x = theCurve->control.x;
  62.         theLine->last.y = theCurve->control.y;
  63.     
  64.     } else if ( (d2 > d1) && (d2 > d3) ) {    // d2 longest
  65.     
  66.         theLine->first.x = theCurve->first.x;
  67.         theLine->first.y = theCurve->first.y;
  68.         theLine->last.x = theCurve->last.x;
  69.         theLine->last.y = theCurve->last.y;
  70.     
  71.     } else {  // d3 must be longest
  72.     
  73.         theLine->first.x = theCurve->control.x;
  74.         theLine->first.y = theCurve->control.y;
  75.         theLine->last.x = theCurve->last.x;
  76.         theLine->last.y = theCurve->last.y;
  77.     
  78.     }//end if
  79.     
  80. }
  81.  
  82. //------------------------------------------------------------------
  83. // Can you say "Pythagore"?
  84. // We convert to floating point for more room to breathe with the
  85. // squares, and for the convenience of sqrt().
  86. double GetLineLength(gxLine *theLine)
  87. {
  88.     double a = (theLine->last.x - theLine->first.x) / kFloatToFixed;
  89.     double b = (theLine->last.y - theLine->first.y) / kFloatToFixed;
  90.     return sqrt(a * a + b * b);
  91. }
  92.  
  93.  
  94. //------------------------------------------------------------------
  95. // "theLength" can be any Fixed number in the range 0x80000000 ... 0x7FFFFFFF
  96. // "theTangent" always is the vector defined by theLine->first ... theLine->last.
  97. void  LineLengthToPoint(double theLength, gxLine *theLine, gxPoint *thePoint, gxPoint *theTangent)
  98. {
  99.     double totallength = GetLineLength(theLine);
  100.     double t;
  101.     Fixed  dx, dy;
  102.     
  103.     if (totallength != 0.0)
  104.         t = theLength / totallength;
  105.     else
  106.         t = 0.0;
  107.     
  108.     dx = theLine->last.x - theLine->first.x;
  109.     dy  = theLine->last.y - theLine->first.y;
  110.     
  111.     if (thePoint) {
  112.         thePoint->x = theLine->first.x + MULTIPLY_DOUBLE_WITH_FIXED(t, dx);
  113.         thePoint->y = theLine->first.y + MULTIPLY_DOUBLE_WITH_FIXED(t, dy);
  114.     }
  115.     
  116.     /* now normalize the tangent */
  117.     
  118.     if (theTangent && (totallength != 0.0)) {
  119.         theTangent->x = FixedDivide( dx, (Fixed)(totallength * kFloatToFixed));
  120.         theTangent->y = FixedDivide( dy, (Fixed)(totallength * kFloatToFixed));
  121.     }//end if
  122. }
  123.  
  124.  
  125. //------------------------------------------------------------------
  126. // The conditions (x(0), y(0)) = theCurve->first, (x(1), y(1)) = theCurve->last,
  127. // (x'(0), y'(0)) = theCurve->control - theCurve->first (tangent at "first")
  128. // (x'(1), y'(1)) = theCurve->last - theCurve->control (tangent at "last")
  129. // are sufficient to determine the coefficients of the quadratic parametrization.
  130. static void ComputeCurveCoefficients(gxCurve *theCurve, CurveCoefficients* cc)
  131. {
  132.     double fixedToFloat = 1 / kFloatToFixed;
  133.     
  134.     cc->a0 = fixedToFloat * theCurve->first.x;
  135.     cc->a1 = fixedToFloat * ((theCurve->control.x - theCurve->first.x) << 1);
  136.     cc->a2 = fixedToFloat * (theCurve->last.x - (theCurve->control.x << 1) + theCurve->first.x);
  137.  
  138.     cc->b0 = fixedToFloat * theCurve->first.y;
  139.     cc->b1 = fixedToFloat * ((theCurve->control.y - theCurve->first.y) << 1);
  140.     cc->b2 = fixedToFloat * (theCurve->last.y - (theCurve->control.y << 1) + theCurve->first.y);
  141.     
  142.     cc->isLinear = ( cc->a2 == 0.0 ) && (cc->b2 == 0.0) ? true : false;
  143. }
  144.  
  145.  
  146. //------------------------------------------------------------------
  147. // That's the integrand sqrt(x'(t)^2 + y'(t)^2) (if you look closely)
  148. static double LengthFunction(CurveCoefficients* cc, double t)
  149. {
  150.     double A = 4 * ((cc->a2 * cc->a2) + (cc->b2 * cc->b2));
  151.     double B = 4 * ((cc->a1 * cc->a2) + (cc->b1 * cc->b2));
  152.     double C = (cc->a1 * cc->a1) + (cc->b1 * cc->b1);
  153.  
  154.     return sqrt(C + t * (B + t * A));
  155. }
  156.  
  157.  
  158. //------------------------------------------------------------------
  159. // After all, this is the most appropriate solution:
  160. // Use the closed formula for the length integral!
  161. // (Mathematica did it for me ...)
  162. static double EvaluateClosedFormula(CurveCoefficients* cc, double x)
  163. {
  164.     double A, B, C, SC, SCBX, result;
  165.     
  166.     A = 4 * ((cc->a2 * cc->a2) + (cc->b2 * cc->b2));    // should always be > 0
  167.     if (A == 0.0)
  168.         {
  169.         DebugStr("\pCurve not quadratic");
  170.         return 0;
  171.         }
  172.     B = 4 * ((cc->a1 * cc->a2) + (cc->b1 * cc->b2)) / A;
  173.     C = ((cc->a1 * cc->a1) + (cc->b1 * cc->b1)) / A;
  174.     SC = sqrt(C);
  175.     SCBX = sqrt(C + x * (B + x));
  176.     result = 4.0 * x * SCBX - 2.0 * B * (SC - SCBX);
  177.     result = result + (B * B - 4.0 * C) * (log(0.5 * B + SC) - log(0.5 * B + x + SCBX));
  178.     result = sqrt(A) * result / 8.0;  
  179.     return result;
  180. }
  181.  
  182.  
  183. //------------------------------------------------------------------
  184. static double NewLengthToParameterFunction(CurveCoefficients* cc, double s)
  185. {
  186.     double a = 0.0;
  187.     double b = 1.0;
  188.     double c, d;
  189.     
  190.     // could set error if s is "out of range"
  191.     
  192.     if (s <= 0.0)
  193.         return 0.0;
  194.     if (s >= EvaluateClosedFormula(cc, 1.0))
  195.         return 1.0;
  196.         
  197.     do    {
  198.         c = 0.5 * (a + b);
  199.         d = s - EvaluateClosedFormula(cc, c);
  200.         if (d < 0)    // c too far right
  201.             b = c;
  202.         else        // c too far left
  203.             a = c;
  204.         }
  205.     while (fabs(d) > kEpsilon);
  206.  
  207.     return c;
  208. }
  209.  
  210.  
  211.  
  212. //------------------------------------------------------------------
  213. double GetCurveLength(gxCurve *theCurve)
  214. {
  215.     CurveCoefficients    cc;
  216.     double    result;
  217.     
  218.     ComputeCurveCoefficients(theCurve, &cc);
  219.     
  220.     if (cc.isLinear) {
  221.     
  222.         gxLine        theLine;
  223.         CurveToLine(theCurve, &theLine);
  224.         result = GetLineLength(&theLine);        
  225.     
  226.     } else {
  227.     
  228.         result = EvaluateClosedFormula(&cc, 1.0);
  229.         
  230.     }//end if
  231.         
  232.     return result;
  233. }
  234.  
  235.  
  236. //------------------------------------------------------------------
  237. static void EvaluateCurveParametrization(CurveCoefficients*    cc, double t, gxPoint* thePoint)
  238. {
  239.     if (thePoint) {
  240.         thePoint->x = (Fixed) (kFloatToFixed * (cc->a0 + t * (cc->a1 + t * (cc->a2))));
  241.         thePoint->y = (Fixed) (kFloatToFixed * (cc->b0 + t * (cc->b1 + t * (cc->b2))));
  242.     }//end if
  243. }
  244.  
  245. //------------------------------------------------------------------
  246. void  CurveLengthToPoint(double theLength, gxCurve *theCurve, gxPoint *thePoint, gxPoint *theTangent)
  247. {
  248.     CurveCoefficients    cc;
  249.     double tx, ty, norm;
  250.     double t;
  251.     
  252.     ComputeCurveCoefficients(theCurve, &cc);
  253.     
  254.     if (cc.isLinear) {
  255.     
  256.         gxLine        theLine;
  257.         CurveToLine(theCurve, &theLine);
  258.         
  259.         LineLengthToPoint(theLength, &theLine, thePoint, theTangent);
  260.         
  261.     } else {
  262.     
  263.         t = NewLengthToParameterFunction(&cc, theLength);
  264.         EvaluateCurveParametrization(&cc, t, thePoint);
  265.         tx = cc.a1 + t * 2.0 * (cc.a2);        // that's the derivative
  266.         ty = cc.b1 + t * 2.0 * (cc.b2);
  267.         norm = sqrt(tx * tx + ty * ty);
  268.         if (theTangent && (norm > 0)) {
  269.             theTangent->x = (Fixed) (kFloatToFixed * tx / norm);
  270.             theTangent->y = (Fixed) (kFloatToFixed * ty / norm);
  271.         }
  272.         
  273.     }//end if
  274. }
  275.  
  276.  
  277. //------------------------------------------------------------------
  278. void DetermineSubLine(gxLine *theLine, double a, double b, gxLine *subLine)
  279. {
  280.  
  281.     LineLengthToPoint(a, theLine, &subLine->first, nil);
  282.     LineLengthToPoint(b, theLine, &subLine->last, nil);
  283.  
  284. }
  285. //------------------------------------------------------------------
  286. void DetermineSubCurve(gxCurve *theCurve, double a, double b, gxCurve *subCurve)
  287. {
  288.     CurveCoefficients    cc;
  289.     double    t0, t1;
  290.     gxPoint    midPoint;
  291.     
  292.     ComputeCurveCoefficients(theCurve, &cc);
  293.     
  294.     if (cc.isLinear) {
  295.     
  296.         /* Treat the curve as a line if it was linear */
  297.         gxLine    theLine, subLine;
  298.         
  299.         CurveToLine(theCurve, &theLine);
  300.                 
  301.         DetermineSubLine(&theLine, a, b, &subLine);
  302.         
  303.         subCurve->first.x = subLine.first.x;
  304.         subCurve->first.y = subLine.first.y;
  305.         subCurve->control.x = (subLine.first.x + subLine.last.x) / 2;  // make it the midpoint
  306.         subCurve->control.y = (subLine.first.y + subLine.last.y) / 2;
  307.         subCurve->last.x = subLine.last.x;
  308.         subCurve->last.y = subLine.last.y;
  309.     
  310.     } else {
  311.     
  312.         /* compute the sub curve */
  313.         
  314.         t0 = NewLengthToParameterFunction(&cc, a);
  315.         t1 = NewLengthToParameterFunction(&cc, b);
  316.         EvaluateCurveParametrization(&cc, t0, &subCurve->first);
  317.         EvaluateCurveParametrization(&cc, t1, &subCurve->last);
  318.         EvaluateCurveParametrization(&cc, 0.5 * (t0 + t1), &midPoint);
  319.         // midPt = (start + end) / 4 + control / 2 ; control = midPt * 2 - (start + end) / 2
  320.         // note that we choose to ignore arithmetic overflow, here: 
  321.         // usable coordinate range limited to -0x3FFF ... 0x3FFF !
  322.         subCurve->control.x = (midPoint.x << 1) - ((subCurve->first.x + subCurve->last.x) >> 1);
  323.         subCurve->control.y = (midPoint.y << 1) - ((subCurve->first.y + subCurve->last.y) >> 1);
  324.         
  325.     }//end if
  326. }
  327.  
  328.